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ABSTRACT 



Signal analysts have traditionally relied on the Discrete 
Fourier Transform and various data windowing schemes for 
signal detection and classification. Some signals, notably 
those of a transient nature, are inherently difficult to 
analyze with these traditional tools. The Discrete Wavelet 
Transform has recently generated considerable interest in 
several areas of digital signal processing and a determination 
of its suitability as a signal analysis tool is necessary. 
Associated with wavelet theory is the concept of 
multiresolutional analyses which allow examination of a signal 
at different scales. 

This thesis investigates dyadic discrete wavelet 
decompositions of signals. A new multiphase wavelet transform 
is proposed and investigated. The • multiphase transform 
technique is shown to be useful in transient signal analysis. 
Several MATLAB programs that perform multiresolutional 
analyses with various supporting features are provided. 
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INTRODUCTION 



The Discrete Fourier Transform has traditionally been the 
dominant tool in digital signal processing for extracting 
waveform features required in the analyses of signals. 
Characteristic of the Fourier transform is the global manner 
in which it operates on the input signal, thereby rendering 
this technique incapable of providing time localization of 
frequency components. The analyses of signals that are 
transitory in nature will particularly suffer the deleterious 
effects of this limitation and the utility of the discrete 
Fourier transform for transient signal analyses may therefore 
be diminished. Various data windows have been used in the 
past in an attempt to improve performance, but all involve 
tradeoffs in resolution and sidelobe levels that may 
unacceptably affect the analysis. [Ref. l:pp. 63-82] 

Recent appearance in the literature of several articles 
describing wavelet transforms and the related multiresolution 
analyses have created interest in determining their 
suitability for signal analysis applications; specifically, to 
overcome the time localization limitation of the Fourier 
transform. This thesis examines the potential of the wavelet 
transform for signal analysis purposes. 
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II. REVIEW OF WAVELET THEORY 



A. DISCRETE WAVELET TRT^SFORM 



Wavelets are families of functions 



ix) 



a,b€ Bt 



( 2 . 1 ) 



a*0 



generated by the dilations and translations of a single 
function. For digital signal processing applications it is 
necessary to discretize the parameter values and fix the 
initial dilation step a^ and the translation step such that 

ao>l, b^*0 

Equation 2.1 becomes. 



The translation parameter is therefore dependent on the 
dilation parameter. A large, positive m will contract the 
function g and cause the translation step to shorten while 
a large negative m dilates the function g and lengthens the 
translation step size. The inverse relationship ensures 
coverage of the entire range. [Ref. 2: pp. 909-910] 



=^o^^h{aoX-nbQ) m,neZ 



( 2 . 2 ) 



with a=ao^ and b=nbQal"' 
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Associated with the discrete wavelet is the discrete 

wavelet transform, T, which maps the function f to a sequence 

indexed by Z^, 

(Tf)^= { h^,f) h{a^x-nb^)fU)dx (2.3) 

If 

[ha)]\^d^ < <-, (2.4) 

and h has sufficient decay, and if T has a bounded inverse on 
its range, the set <h^> forms a frame for all square- 
integrable, one-dimensional functions f, 

fix)€L^{?lt) (2.5) 

The significance of a frame is that numerically stable 
algorithms exist which allow f to be reconstructed from the 
wavelet coefficients <h^,f>. [Ref. 2: p. 911] 

Selection of appropriate h,ag, and bg is influenced by 
various factors. The desire to take advantage of previously 
published work and to minimize redundancy in the wavelet 
representation resulted in choices of ag = 2, bg = 1, and 
several h (discussed below) such that the constitute an 
orthonormal basis of compact support. 
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B. MULTIRESOLUTION ANALYSIS 



Multiresolution analysis, as the name implies, describes 
the L function / as a series of approximations of the 
function at different levels of resolution and consists of a 
family of embedded closed subspaces 

(Si) : ... cV.2cV_iCV’oCV;^cV2 ... (2.6) 

such that 
and 

f(-) 6V„ ~ (2.8) 

Also, there is a scaling function <|)(x)eVo such that the 
(|)^ (where (j)^(x) =2'"^^<|) (2'”x-n) ) constitute a basis for 

= span((t>^) (2.9) 

Let P^f be the orthogonal projection of f onto V^. From the 
preceding, it can be seen that liin^_„P^f = f, for all f€L^{R) . 

[Ref. 3:pp. 967-968] 
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Considering our parameter choices and the decision to use 
orthonormal bases, if we define 

( 2 . 10 ) 

then 

V = ^ (2.11) 

Now define Q^f = and as the orthogonal complement 

of so that = Then Ojj, is the 

orthogonal projection of any f e (R) onto the 

orthogonal complement of in . The IV^ are scaled 

versions of , 

f(‘) ehr^^f(2‘”’) e Wq (2.12) 

. 2 

and the are orthogonal spaces which sum to I<“(R) 

l2(R) = © (2.13) 

Similar to the V^, there exists in a wavelet function ^ 
such that 
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= spaFTO 



(2.14) 



where i(r^(x) = (2'"x-n) . If we define 

^mn (f) = f) (2.15) 



then 

Qf = Y'd if) \i! EW (2.16) 

Naturally, the function may be reconstructed from its 
projections onto the bases vector spaces. [Ref. 2:pp. 916-918] 
1. Properties of the Wavelet-Scaling Function Pair 

Following is a brief summary of the relationship 

between the (J) and i|f families. Since <J)q (x) eVoCi/^ , it can 
be written as 



ee 

4>o(^) = )) (2.17) 

n—~ ^ ^ 
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Let 



h{n) = J 4>o (x)<l)o (2 x- 73) dx 

= 2"^(<l)o(x) ,4>,(x-|)) 



( 2 . 18 ) 



so that 

oo 

<t>o(x) = ^/2Y, h(n)<t>,(x-^) 

OQ 

= 2 h(jn)<J)o(2j!c-j2) 



The Fourier transform of this equation is 

(2q) = //(co) O (co) 



( 2 . 19 ) 



( 2 . 20 ) 



where we have defined 

H(o>) = A (n) (2.21) 

n=-oo 



To satisfy the above relationships, the following properties 
must hold: 



IH(0)\ = 1 

|H(co)|2 + ih(o+u)|2 = 1 



( 2 . 22 ) 
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From Equation 2.20 it can be seen that the 



4> can thus be 



derived if H(o)) is known since 



p=i 


(2.23) 


Knowledge of (j) can subsequently be used to find 
(x) ^ Wq c , it can be written as 


. Since 


oo 

21=-» ^ ^ 


(2.24) 



Let 

oo 

gin) = (x) 4>o ( 2x-n) 

— oo 
-1 

= 2 2(i|fg(x) ,4>^(x--^)) 


(2.25) 



and by a similar approach to that used previously we can find 



i|»o(x) = g(n)^^ix-^) 

oo 

= 2 g(n) 4>o (2x-n) 

j3 = -oo 


(2.26) 
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The Fourier transform of Equation 2.26 is 



T (2o) = G(o) <I> (o>) 



(2.27) 



The following properties result: 



1G(0)| = 0 

|G(o)|2 + |G(o+n)|2 = 1 



(2.28) 



G((j>) + H(,<ji+n) G(o+7i) = 0 



Equations 2.22 and 2.28 are recognized as characteristic of 
'•conjugate" filters and are necessary to ensure orthogonality 

among the 4> and ijr families of functions. Equation 2.28 
ensures that each function of the 4> (commonly referred to 

as the scaling function) is orthogonal to each function of the 

family. This last condition results from the fact that 
V^± . [Ref. 4: pp. 16-23] 

An example of a function that fulfills the stated 
conditions is 



G(o) = 



(2.29) 
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from which we can find 



gin) = i-l)^-"hil-n) (2.30) 

Equation 2.30 is the defining relationship for describing 
wavelet-scaling function dependency in the application 
algorithms. [Ref. 3: p. 679] 

2 . Decomposition and Reconstruction Algorithms 

Recalling that =(4>j^,f) , we can derive the 

decomposition recursion algorithm: 

= f<t>^.i(x-2-^^-^^n)f(x)dx 

= f h(k) 4>jji(x-2'^^'^^n-2~^k) ] f (x) dx 

* (2.31) 

= v/2j]h(ic) U^{x-2-’”i2n+k))f(x)dx 

k •' 

= s/T^hik) c^{2n+k) 

k 

Likewise we find 

= sf2Y,9ik) cj2n+k) (2.32) 

k 
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To determine the reconstruction recursion algorithm recall 



OP OP 

~ ^ (m-l) k 



k»-“ 



k»-» 



(2.33) 

(2.34) 



and substitute terms to get 



~ ^ ^m-1 ^m-l ( ^4>mn > ^ (m-1) J(} 

k k 



(2.35) 



Using a change of variables the following can be found 





= s[2h (n-2k) 


(2.36) 




= sf2g{n-2k) 


(2.37) 



so that 

c^(n) = / 2 [J^c^.,(k)J 2 (n- 2 Jc) +'£d^.,(Jc)g(n-2Jc)] (2.38) 

Jc k 



The equations of this section demonstrate the pyramidal 
structure of multiresolution analysis and readily lend 
themselves to digital signal processing applications. [Ref 4: 
pp. 25-28] 

The important conclusion is that, assuming knowledge 
of the g and h vectors, the c and d coefficients at any level 
can be completely determined from the c coefficients at the 
next higher level; also, the c coefficients at any level can 
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be determined from the c and d coefficients at the adjacent 
lower level. Notice that calculations within this pyramidal 

algorithm do not require explicit use of the 4> and ijr 

functions since interlevel relationships are defined entirely 
by the g and h vectors. Obviously, the same results would be 
reached if dilations and translations of the orthogonal basis 
functions were used directly at each level in the 
decomposition and reconstruction, but would be found at great 
computational cost because of the requirement to calculate 
inner products. 
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III. COMPUTER IMPLEMENTATION OF MULTIRESOLUTION ANALYSIS 
A. INTRODUCTION 

Our goal is to investigate the use of multiresolution 
analyses to extract information from an input signal. In this 
case, a one-dimensional data sequence will represent the 
input. To avoid having to numerically evaluate any inner 
products we will equate the Cp(n) coefficients with the data 

sequence, thereby constructing an auxiliary function f, with 
/=2^C(j (n) <t)o (n) , which clearly resides in V^q. Since CQ(n) 
entirely describes the input, it represents the highest 

resolution level possible and the apex of the pyramidal 
algorithm. The multiresolution framework previously described 
can now be used to decompose f (thus the data sequence C(j(n)) 
into lower resolution, i.e. m < 0, approximation coefficients 
c^(n) , and detail coefficients d^(n) . The reconstruction 
recursive equations can be used to recover the original f (and 
therefore the data sequence) . A different multiresolution 
analysis exists for each scaling function/wavelet set. This 
thesis is limited to the Haar wavelet and to Daubechies' group 
of compactly supported wavelets. The h vectors for these 
wavelets were published in Reference 1 and are listed in 
Appendix H for convenience. Describing the h vectors is the 
common method of defining scaling function/wavelet sets 
because all other desired information can subsequently be 
derived from them. 
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Several computer programs were written to implement 
various multiresolution analyses. The programs were written 
in MATLAB to take advantage of the transportability of m- 
files between various operating systems. The capabilities of 
these programs, listed in Appendices B-H, will be described in 
general. Specific information can be obtained by referring to 
program comments in the appropriate appendix. 

The normal entry point into the collection of programs is 
the calling program wavelet. m (Appendix B) . A brief 
description of the various possible options is presented upon 
entry and the user can make a decision based on the choices 
provided. Specifically, the user may decide to use the Haar 
wavelet, his own wavelet, or one from the Daubechies group of 
nine compactly supported wavelets and analyze with either a 
single phase or multiphase approach. 

The single phase approach implies the decomposition start 
point is the data start point, i.e. a "snapshot" of the entire 
data sequence is processed. The term multiphase refers to 
starting a decomposition at every possible data 
sequence/wavelet phase relationship, i.e. a "sliding window" 
approach, to maximize the signal analysis capability. As 
desired, wavelets are phase sensitive but feature extraction 
may be hindered if the "best" input phase relationship is not 
used. The multiphase analysis therefore contains multiple 
sets of coefficients, with each set containing all the 
information necessary for reconstruction. 
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B. SINGLE PHASE MULTIRESOLUTION ANALYSIS 

The second program, haarwave.m (Appendix C) , does the 
"snapshot" analysis of an input data sequence using the Haar 
basis. Because of the simple nature of the Haar, it is easy 
to generate the actual approximation and detail function at 
each resolution level. The input data is viewed as a 
piecewise constant function (equating to the sample and hold 
operation) to allow easy inner product calculation. 
Representations of the properly weighted scaling function and 
wavelet at each level are provided to clearly indicate the 
dyadic relationship between levels. Additionally, the 
reconstruction of the original sequence can be viewed as well 
as a plot of the reconstruction error. Finally, distribution 
of the energy of each approximation and detail coefficient and 
the total energy of each resolution level are displayed. 
These energy distributions are the foundation for signal 
analysis. 

The third program, daubwave.m (Appendix D) , allows the 
user to pick one of the wavelets from Daubechies group or to 
enter a valid user-defined h vector to define the basis set. 
Basically, the same output plots that were generated for 
wvhaar.m can be seen. Because of the complexity and 
irregularity of the basis functions, however, only the value 
of the coefficients which are generated for specific points 
will be plotted instead of the inner product representation. 
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C. MULTIPHASE WAVELET ANALYSIS 

The last analytical program, multiphs.m (Appendix E) , 
allows the use of any of the previously defined wavelets with 
the multiphase algorithms. The different coefficient sets are 
displayed simultaneously, which is possible because their 
coefficient indices interleave without interference. Clearly, 
redundant information is provided in the sense that there are 
more coefficients than would be necessary for reconstruction 
of the data sequence over its interval of support. However, 
because of the different zero padding necessary for each set 
of coefficients a different P^f in is defined for each set. 
The user will therefore be able to see the P^f with the most 
distinctive set of coefficients to ease analysis. 

D. SUPPORTING PROGRAMS 

The program basisplt.m (Appendix F) supports the 
daubwave.m and multiphs.m programs by generating iterative 
plots of the scaling function and wavelet bases. It can also 
be called up directly and will use the vector "hcoeff as the 

h coefficients to determine the 4> and i|f based on a 
graphical recursion method. 

The functions wvinput.m (Appendix G) and daubdata.m 
(Appendix H) are called as necessary by any of the 
multiresolution analysis programs to provide selected input 
data signals or the h vectors from Daubechies' group, 
respectively. The input signal options consist of a sinusoid. 
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a pseudo-noise (PN) sequence, a sinusoid modulated with a PN 



sequence, 

parameters 



and any of the above corrupted by noise. Various 
can be modified to satisfy the user's needs. 
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IV. SIGNAL ANALYSIS 



All plots associated with this section can be found in 
Appendix A. 

A. SINGLE PHASE ANALYSIS OF A SINUSOIDAL WAVEFORM 

The first signal to be analyzed will be the sinusoid shown 
in Figure 1, along with its level 0 approximation. Figures 2- 
5 show the Haar multiresolution analysis, stopping at level -4 
since there is no energy at any lower level for this 
decomposition. The reconstruction plots (not shown) are 

indistinguishable from the decomposition plots because the 
maximum reconstruction error is 15 orders of magnitude below 
signal levels. The energy contained at each resolution level 
for both approximation and detail coefficients can be seen in 
Figure 6. Clearly, the maximum energy change occurs in the 
detail coefficients at level -4 and the energy in the 
approximation coefficients at every level is the sum of the 
energy of the detail and approximation coefficients from the 
level below, as expected. Figures 7 and 8 provide a clear 
summary of exactly where the signal and filter best match with 
respect to sample number and resolution level . In this 
example, the Haar decomposition highlights the phase changes 
associated with the switching between positive and negative 
half-cycles because of the signal/wavelet phase coherence. 
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The decomposition is repeated with a different 
multiresolution analysis based on the 6-coefficient (third 
order) Daubechies wavelet. Only the coefficient energy plots 
are shown in Figures 9-11 for comparison with the Haar plots. 
It can be shown that the time index of coefficients goes 
beyond the original support length of the data (potentially 
out to sample 257 in this example, although the plots were 
truncated at sample 90 for display purposes and because very 
little energy is associated with samples beyond 90) . Although 
the coefficients outside the signal support range are 
generally very small, they are necessary for completeness. 
Computing these coefficients requires affixing zeros to the 
sequence, most significantly at the trailing edge. These 
"edge effects” are minimized at the leading edge (no 
coefficients are generated for negative time indices) by 
shifting the h vector so that nonzero values may initially 
exist only over [- (m-2) , . . . , 0 , 1] and then assigning the 
coefficient value to the index of the second term from the 
right (initially 0). This is intuitively satisfying from a 
causality viewpoint and also aids interpretation by only 
having coefficients outside the data sequence on one side. 
Leading edge effects consist of the influence of leading zeros 
necessary for filter coefficients with negative time indices. 
Note that edge effects increase as the resolution level 
decreases because lower resolution coefficients are influenced 
by a larger number of data points. The extent of edge effects 



19 



is also dictated by the number of h coefficients. 
Coefficients beyond the data sequence support will not have 
time indices greater than [ (2 "’-1) (number of h coefficients - 
1) ] for any of the decomposition programs and often 
significantly fewer (e.g., if the data sequence length is a 
power of 2 the Haar decomposition will not generate any terms 
beyond the data length, as can be seen in the previous plots) . 

The next set of plots (Figs. 12-18) show the phase 
sensitivity of the decomposition when the input sinusoid phase 
is shifted from 0 to 45 degrees. The energy distribution is 
no longer as concentrated as in the previous set of plots and 
interpretation is consequently more difficult. Changing the 
sampling frequency relative to the sinusoid frequency will 
have a similar effect. 

B. SINGLE PHASE ANALYSIS OF A BINARY PHASE SHIFT KEYED SIGNAL 

The Binary Phase Shift Keyed signal shown in Figure 19 
will serve to highlight a deficiency in the single phase 
approach. Notice with the Haar wavelet in Figures 20-22 that 
the decomposition is exactly the same as the pure sinusoid, 
i.e., because of the phase alignment we cannot discriminate 
between the two signals' coefficient energy plots. The 
Daubechies 6-coefficient wavelet decomposition (Figs. 23-25) 
does appear significantly different than the sinusoid 
decomposition but there is no clear indication of every phase 
reversal. Thus, the results are ambiguous. 



20 



C. MULTIPHASE ANALYSIS OF A BINARY PHASE SHIFT KEYED SIGNAL 



Since these responses would be unsuitable for signal 
analysis purposes the multiphase approach (described 
previously) was developed. It can be viewed as starting 
decompositions at each of the 2’"' phases of each level so that 
the most distinctive representation can be used for analysis 
regardless of initial signal phase. The multiphase plots are 
shown in Figures 26 and 27 for the Haar and in Figures 28 and 
29 for the Daubechies 6-coefficient filter. There is clear 
indication in both sets of plots of the signal's phase 
inversions . 

D. MULTIPHASE ANALYSIS OF A TRANSIENT SIGNAL 

Finally, the transient signal of Figure 30 was 
investigated. Use of the "zoom" feature in the programs was 
necessary to clearly see the response. Figures 31-34 show the 
responses associated with the Haar and Daubechies 6- 
coefficient wavelets that have been the standard throughout 
this thesis. Higher order Daubechies (more regular) wavelets 
were also experimented with in an attempt to best fit the 
signal and it can be seen that the 18-coefficient (ninth 
order) wavelet used for Figures 35 and 36 does a good job of 
clearly indicating the presence of the signal and of 
maximizing the concentration of energy into only a few points 
in a single resolution level. This good "match" of the signal 
with the wavelet filter could be used for signal 
classification purposes as well as detection if an a priori 
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selection of the wavelet filter based on a similar known 
signal had occurred. 
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V. CONCLUSIONS 



Wavelet-based loultiresolution analysis is another tool for 
the signal analyst and great potential exists for its use in 
various applications. Signals that can only be represented in 
the frequency domain with large numbers of significant terms 
provide special motivation for wavelet-based decomposition, as 
can be seen in the transient signal example above. Signal 
detection and pattern recognition through the use of "matched" 
wavelets may also prove useful, especially in areas where 
analysis at different resolution levels (i.e., possible signal 
dilation/contraction) is required. The extent of the ultimate 
usefulness and popularity of wavelets will become known as 
others gain knowledge of basic wavelet characteristics and 
incorporate multiresolution analyses into their research. 
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APPENDIX A 



SIGNAL ANALYSIS PLOTS 
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MATLAB Plot of Input Signal Vector 




Figure 1. 



Input Sinusoid and Approximation 
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Sample number "n" (Time = n * 0.0625 sec) 



Approximation function at resolution level - 




Figure 2. Haar Response at Resolution Level -1 
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Sample number ”n” (Time = n * 0.0625 sec) 



Approximation function at resolution level -2 
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Figure 3. Haar Response at Resolution Level -2 



27 



Sample number "n" (Time = n * 0.0625 sec) 



Approximation function at resolution level -3 
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Figure 4. Haar Response at Resolution Level -3 
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Sample number "n" (Time = n * 0.0625 sec) 



Approximation function at resolution level -4 
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Figure 5. Haar Response at Resolution Level -4 
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Sample number ’*n" (Time = n * 0.0625 sec) 



Normalized Energy of Approximation Function vs. Resolution Level 





Figure 6. 



Haar Resolution Level 



Energy Distribution 
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Resolution Level 



Energy in Approximation Coefficients using Haar basis 




Figure 7. Haar Approximation Coefficient Distribution 
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Energy in Difference Coefficients using Haar basis 




Figure 8. 



Haar Detail Coefficient 



Distribution 
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Normalized Energy of Approximation Coefficients 





Figure 9. Daubechies Resolution Level Energy Distribution 
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Resolution Level 



Energy in Approximation Coefficients 




figure 10. Daubechies Approximation CoeHlcient Distribution 
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ine Daubechies Wavelet of order 3 



Energy in Detail Coefficients 




Figure 11. 



Daubechies Detail 



Coefficient Distribution 
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ing Daubechies Wavelet of order 3 



MATLAB Plot of Input Signal Vector 
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Sample number *'n" (Time = n * 0.0625 sec) 



Normalized Energy of Approximation Function vs. Resolution Level 





Figure 13. Haar Resolution Level Energy Distribution 
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Resolution Level 



Energy in Approximation Coefficients using Haar basis 




Figure 14 . Haar Approximation Coefficient Distribution 
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Energy in Difference Coefficients using Haar basis 




Figure 15. 



Haar Detail Coefficient Distribution 
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Normalized Energy of Approximation Coefficients 
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Figure 16, Daubechies Resolution Level Energy Distribution 
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Resolution Level 



Energy in Approximation Coefficients 




Figure 17. Daubechies Approximation Coefficient Distribution 
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using Daubechies Wavelet of order 3 



Energy in Detail Coefficients 




Figure 18. 



Daubechies Detail 



Coefficient Distribution 
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ing Daubechies Wavelet of order 3 



MATLAB Plot of Input Signal Vector 




Figure 19. Binary Phase Shift Keyed Signal and Approximation 
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Sample number ”n" (Time = n * 0.0625 sec) 



Normalized Energy of Approximation Function vs. Resolution Level 





Figure 20. Haar Resolution Level Energy Distribution 
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Resolution Level 



Energy in Approximation Coefficients using Haar basis 
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Figure 21. 



Haar Approximation Coefficient Distribution 
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Energy in Difference Coefficients using Haar basis 




Figure 22. 



Haar Detail Coefficient 



Distribution 
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Normalized Energy of Approximation Coefficients 





Figure 23. Daubechies Resolution Level Energy Distribution 
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Resolution Level 



Energy in Approximation Coefficients 




Figure 24. 



Daubechies 



Approximation Coefficient Distribution 
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using Daubechies Wavelet of order 3 



Energy in Detail Coefficients 




Figure 25. 



Daubechies Detail 



Coefficient Distribution 
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inc Daubechies Wavelet of order 3 



Energy in Approximation Coefficients 
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Figure 26 . 



Multiphase Haar 



Approximation Coefficients 
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Energy in Oetnil Coefficients 




Figure 27 . 



Multiphase 



Haar Detail Coefficients 
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Energy in Approximation Coefficients 




Figure 28. Multiphase Daubechies Approximation Coefficients 
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using Daubechies Wavelet of order 3 



Energy in Detail Coefficients 




Figure 29. 



Multiphase Daubechies Detail Coefficients 



53 



using Daubechies Wavelet of order 3 



MATLAB Plot of Input Signal Vector 




Figure 30. Transient Signal 
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Sample number 



Energy in Approximation Coefficients 




Figure 31. Multiphase Haar Approximation Coefficients 
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Energy in Detail Coefficients 




Figure 32. Multiphase Haar Detail Coefficients 



56 



Energy in Approximation Coefficients 




Figure 33. Multiphase Daubechies Approximation Coefficients 
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using Daubechies Wavelet of order 3 



Energy in Detail Coefficients 




Figure 34. 



Multiphase Daubechies Detail Coefficients 
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lies Wavelet of order 3 



Energy in Approximation Coefficients 




Figure 35. Multiphase Daubechies (Order 9) Approx. Coefficients 



59 



using Daubechies Wavelet of order 9 



Energy in Detnil Coefficients 




Figure 36. 



Multiphase Daubechies 



(Order 9) Detail Coefficients 
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inc Daubechies Wavelet of order 9 



APPENDIX B 



CALLING PROGRAM FOR ACCESSING MULTIRESOLUTION ANALYSES 

PROGRAMS AMD FUNCTIONS 
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% Calling program for wavelet multiresolution analysis. The capabilities 
% and options of the entire set of programs are presented, in general terms, 
% in the comments below. 



b=[' 
N1 = 






N2 = [ ' 



This program will perform multiresolution analysis on input data' 
through the use of various orthonormal bases of compact support. ' 

It is interactive in nature and will allow you the opportunity ' 
to select among several different options to best support your ' 
specific needs. Use the Return key to control movement through ' 

the program (after prompts and questions, to view next plot, etc.)' 

/ 

Program options include: ' 

/ 

1. Selection of Basis Functions ' 

a. Haar wavelet ' 

b. Daubechies group of compactly supported wavelets ' 

c. User-input wavelet ("h" coefficients) ' 

/ 

2. Generation of Scaling Function/Wavelet plots ' 

t 

3. Selection of Input Data ' 

a. User-defined data vector (with optional noise background): ' 

1) Sine wave ' 

2) Pseudo-noise (PN) sequence ' 

3) Sine wave modulated by PN sequence ' 

b. User-input data vector ' 

<Return> for more ...']; 

4. Decomposition algorithm based on: ' 

a. Data "snapshot" - Fixed start point; one decomposition; ' 

single phase approach ' 

b. "Sliding" Data window - Floating start point; decomposition ' 

for each start point; multiphase ' 

approach; discrete approximation to ' 

continuous wavelet transform ' 



Do you want to see: ' 

/ 

1. Haar wavelet data "snapshot" decomposition, ' 

(input data treated as a piecewise constant function, ' 

allowing inner products to be calculated and displayed) ' 

/ 

2. Daubechies or user wavelet data "snapshot" decomposition, or ' 

/ 

3. "Sliding" data window wavelet (any type) decomposition? '] 



clc 

disp (b) 
disp(Nl) 
pause 
clc 

disp (b) 
disp (N2) 
disp (b) 

desire = input ( 'Answer 1, 2, or 3 . '); 
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1 



if desire == 
haarwave 
elseif desire 
daubwave 
else 

multiphs 

end 



APPENDIX C 



PROGRAM FOR SINGLE PHASE ANALYSIS USING HAAR WAVELET 
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% This program performs single-phase Haar wavelet decomposition on a 
% user-supplied or program-generated input waveform. Representation is 
% made with bar graphs to show actual inner product calculation (input 
% data is taken as a piecewise constant function) . The input data sequence 
% will be zero-padded to accomodate the calculation of all possible non- 
% zero coefficients and approximation ("c”) and detail ("d*') coefficients 
% may be computed for sample numbers beyond the data endpoint. The user 
% should consider these "edge effects" in his analysis. 

% 

% Determine if user wants to input an external data vector or desires to 
% build one through the program. 

clc 



b=[' 

Q1 = ['Do you wish to use: ' 

' 1. your own MATLAB formatted row vector, or ' 

' 2. a program generated vector from the following menu?' 

' - Sine wave ' 

' - Pseudo-noise (PN) sequence ' 

' - Sine wave modulated by PN sequence 

disp(b) 

disp(Ql) 

Pick = input ( 'Answer 1 or 2 : '); 

if Pick == 1, 

% Read in user's input vector 

N1 = [ ' Note; If the number of data points is a power of 2, the' 

' results are easier to interpret because there are not ' 

' any "edge effects". ']; 

disp(b) 
disp(Nl) 

Wavedata = input ('What is the name of your input vector? '); 
sampfreq = input ('What v/as the sampling frequency (Hz)? '); 

Tsample = 1/sampfreq; 

else 

[Wavedata , Tsample] = wvinput (Pick) ; 

end 



% Decomposition algorithm h(0)=0.5, h(l)=0.5, g(0)=0.5, g(l)=-0.5 
datlngth - length (Wavedata) ; 

numrows = ceil (log(datlngth)/log(2) ) ; % find number of resolution levels 

numpts = 2^numrows; Newdata = zeros ( 1 , numpts) ; 

Newdata ( 1; datlngth) = Wavedata; 

d = zeros (numrows , numpts) ; c = zeros (numrows+1 , numpts) ; 
detail = zeros (numrows , numpts) ; approx = zeros (numrows , numpts) ; 
c ( numrows+ 1 , : ) = Newdata ; 
top = numpts ; 

const = l/sqrt(2) ; % normalization constant * coefficient magnitude 

ctr = numrows ; 
while ctr >= 1, 

n = numrows-ctr+1 ; 
shift = 2^n; 
shiftl = shift/2; 

top = ceil (top/2); % number of points at this level 
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% get "c" and "d" coefficients 



end 



for k = l:top 

jump =shift*k; 
indexl = jump- (shift-1) ; 
index2 = jump- (shiftl-1) ; 
firsttrm = c (ctr+1 , indexl) ; 
if index2 <= numpts, 

secndtrm = c (ctr+1 , index2) ; 
else 

secndtrm == 0 ; % zero padding 



end 

d (ctr, indexl) = firsttrm 
c(ctr, indexl) = firsttrm 
for j = 0; (shift-1) 
if j < shiftl, 

detail (ctr , indexl+j ) 
else 

detail (ctr, indexl+j ) 
end 

approx (ctr, indexl+j ) = c (ctr, indexl) ; 
end 
end 

d(ctr,:) = const*d (ctr , : ) ; 
c(ctr,:) = const*c (ctr , : ) ; 
normlize = sgrt (2*shift ) ; 
detail (ctr ,: ) = detail (ctr ,:) /normlize ; 
approx (ctr ,: ) = approx (ctr ,:) /normlize ; 
ctr = ctr-1; 



secndtrm; 

secndtrm; 

% build matrices for display purposes 
d (ctr , indexl) ; 

-d (ctr , indexl) ; 



% normalization constant for this level 



% Plot "c" and "d" coefficients by resolution level 
indxtoO = [ 0 . 5 : numpts-0 . 5 ] ; 

plotmin = 1 . 2*min (Wavedata) ; plotmax = 1 . 2 *max (Wavedata) ; 

if ( (plotmin==0) & (plotmax==0) ) , plotmax = 0.5; end 

if plotmin > 0, plotmin = 0; end 

if plotmax < 0, plotmax = 0; end 

v= [ 0 , numpts, plotmin , plotmax] ; 

axis (V) ; 

bar ( indxtoO , c (numrows+1 , : ) ) 

title ( 'Approximation function to Input Signal at resolution level O') 
xlabel ([ 'Sample number "n" (Time = n * ' , num2str (Tsample) , ' sec)']) 
pause 
clg 

outptdet = zeros (1,1); outptapp == zeros (1,1); 
for k = numrows:-!:! 

level = k-numrows-1; 
step = (2^ (-level) )/2 ; 
for j = 1 : numpts/step 

outptdet(j) = detail (k, (j-1) *step+l) ; 

end 

stepa = 2*step; 

indxtoOa = [ stepa/2 : stepa : numpts-stepa/2 ] ; 
for j = 1 : numpts/stepa 

outptapp(j) = approx (k, (j-1) *stepa+l) ; 

end 

axis (v) ; 
if k > 1, 

subplot (211) , bar (indxtoOa, outptapp) 
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else 

indxtoOa = [0:numpts]; 

approxl = approx (1, :) ;approxl (nuiupts+1) = approx(l,numpts); 
subplot (211) , plot (indxtoOa, approxl) 

end 

title ([ 'Approximation function at resolution level num2str ( level )] ) 
xlabel ([ 'Sample number *'n” (Number of samples = ' , num2str (numpts) , ' ) ' ] ) 
axis (v) ; 

subplot (212) , bar ( indxtoO , outptdet ) 

title ([' Detail function at resolution level ', num2str (level )} ) 
xlabel ([ 'Sample number "n*’ (Time = n * ' , num2str (Tsample) , ' sec)']) 
pause 
clg 

clear outptapp outptdet indxtoO 
indxtoO = indxtoOa; 
clear indxtoOa 

end 

subplot 

% 

% Reconstruction algorithm 

clc 

recmpout = zeros (1,1); 
disp(b) 

N2 = [ ' Level-by-level Recomposition can be observed if desired. ' 

' The Recomposition algorithm starts with the lowest level ' 

' Approximation Function and successively adds in the Detail' 

' Function to obtain the next higher level Approximation. ']; 
disp (N2) 
disp (b) 

ckalgthm = input ('Do you want to see the Recomposition (Y or N) ? ','s'); 
if ckalgthm == 'Y', 
level = -numrows; 
axis (v) ; 

plot (indxtoO , approxl) 

title ([' Level ', num2str (level ), ' Recomposition']) 

xlabel (' Sample number') 

pause 

clg 

recomp = approx (1, :) ; 
for k = 1: numrows 
level = level+1; 
recomp = recomp+detail (k, : ) ; 
step = 2^ (-level); 

indxtoO = [step/2 : step: numpts-step/2] ; 
for j = 1 : numpts/step; 

recmpout (j) = recomp ( (j -1) *step+l) ; 

end 

axis (v) ; 

bar ( indxtoO , recmpout) 

title ([ 'Level ', num2str (level) , ' Recomposition']) 

xlabel ( 'Sample number') 

pause 

clg 

end 

bar ( indxtoO , Newdata-recmpout) 
title ( 'Recomposition Error') 
xlabel (' Sample number') 
pause 
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clg 

end 



% Coefficient energies are used as a measure of response 

c = c . ^ 2 ; 
d = d.^2; 

Enrgytot = sum (c (numrows+l , ; ) ) ; 

Enrgycro = zeros ( 1 , numrows+1) ; Enrgydro = zeros ( 1 , numrows+l) ; 

for j = l:numrows % find energy in each resolution level 

Enrgycro(j) = (sum(c (j ,:))) /Enrgytot ; 

Enrgydro(j) = (sum (d (j ,:))) /Enrgytot ; 

end 

if Enrgytot == 0, 

Enrgycro (numrows+i) = o; 
else 

Enrgycro (numrows+1) = 1; 
end 

Enrgydro (numrows+1) = 0; 

Ivl = [ - (numrows) : 1 : 0] ; 

V = [- (numrows+1 ), 1 , 0 , 1 . 2 ] ; 
axis (V) ; 

subplot(211) , bar (Ivl , Enrgycro) 

title ( 'Normalized Energy of Approximation Function vs. Resolution Level') 
xlabel ( 'Resolution Level') 
axis (V) ; 

subplot(212) , bar (Ivl , Enrgydro) 

title ( 'Normalized Energy of Detail Function vs. Resolution level') 

xlabel ( 'Resolution Level') 

pause 

subplot 

clc 

% 

% Display mesh and contour plots of coefficient energy with optional "zoom” 

% capability to aid analysis 

strtsarop = 1; endsamp = numpts; strtrow = 1; endrow = numrows; zoom = 1; 
highres = -1; lowres = -numrows; 
while zoom == 1, 

rangea = [ -highres-1 : -lowres] ; ranged = [ -highres : -lowres] ; 

indxtoO = [ strtsamp-1 : endsamp-1] ; 

mesh (c (strtrow: endrow+1 , strtsamp : endsamp) ) 

title ( 'Energy in Approximation Coefficients using Haar basis') 
pause 

contour (c (strtrow: endrow+1 , strtsamp : endsamp) , 10, indxtoO , rangea) 
title (' Contour Map of Approximation Coefficient Energy Distribution') 
xlabel ([' Sample number "n" (Time == n * ' , num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= -Resolution Level)') 
pause 

mesh (d (strtrow: endrow, strtsamp: endsamp) ) 

title (' Energy in Difference Coefficients using Haar basis') 
pause 

contour (d (strtrow: endrow, strtsamp : endsamp) , 10, indxtoO , ranged) 
title (' Contour Map of Difference Coefficient Energy Distribution') 
xlabel ([ 'Sample number "n" (Time = n * ', num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= -Resolution Level)') 
pause 
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clc 

shovmore ~ input('Would you like to "zoom" in on a section (Y or N)? 
if showmore == 'Y', 
disp (b) 

disp(' You will set the display sample number and resolution level limits.') 
disp (b) 

disp([' Sample range 0 : ' , num2str (numpts-1) ] ) 

disp([' Resolution range -1 : ' , num2str (-numrows) ] ) 

revu = input ('Need to see the original contour plots again (Y or N)? ','s'); 
if revu == ' Y' , 

contour ( c , 1 0 , [ 0 : numpts-1 ] , [ 0 : numrows ] ) 

title (' Contour Map of Approximation Coefficient Energy Distribution') 
xlabel ([' Sample number *'n” (Time - n * ' , num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= -Resolution Level)') 
pause 

contour (d, 10, [0:numpts-l] , [ 1 : numrows] ) 

title (' Contour Map of Difference Coefficient Energy Distribution') 
xlabel ([' Sample number *'n" (Time = n * ' , num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= -Resolution Level)') 
pause 
end 

strtsamp = input ( 'Sample start point? ')+!; 
endsamp = input (' Sample endpoint? ')+!; 

highres = input ( 'Highest resolution level (least negative)? '); 
lowres = input (' Lowest resolution level (most negative)? '); 
endrow = numrows+l+highres ; 
strtrow - numrows+l+lowres ; 
else 

zoom = 0; 
end 
clg 
end 
clc 
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APPENDIX D 



PROGRAM FOR SINGLE PHASE ANALYSIS USING 
DAUBECHIES OR USER WAVELET 
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% This program performs single phase decomposition of the input waveform 
% with either a Daubechies wavelet of user-defined order or with a wavelet 
% provided by the user. The input data sequence will be zero-padded to 
% accomodate the calculation of all possible non-zero coefficients and 
% approximation (*'c") and detail ("d") coefficients may be computed for 
% sample numbers beyond the data endpoint. The user should consider these 
% "edge effects" in his analysis. 



% 

% Input the "h" coefficients 

clc 

b = [ ' ' ] ; 
disp(b) 

D1 = [ ' Would you like to use: ' 

' 1. your own decomposition coefficients (use wvhaar.m if you desire' 

' to use the Haar basis set) ' 

' 2. those associated with Daubechies compactly supported wavelets? ']; 

disp(Dl) 

hpick = input ( 'Answer 1 or 2 : ' ) ; 

disp (b) 

if hpick == 1, 

disp(' Input as a row vector with sum of coefficients normalized to "1".') 
hcoeff = input ('What is the name of your decomposition coefficient vector? '); 
choice = 1; 
else 

choice = input ('What is the desired wavelet order (2-10 are available)? '); 
[hcoeff] = daubdata (choice) ; 
end 

hlength = length (hcoeff ) ; 

if hlength == 0 
disp (b) 

disp('ERROR: Wavelet order selected is outside allowable range.') 
break 
end 

% 

% Determine if user wants to see plots of basis functions 

ckplt = input ('Do you want to see the Scaling Function/Wavelet (Y or N)? ','s'); 
if ckplt == 'Y' 
basisplt 
end 

% 

% Get the input data vector 

clc 

disp (b) 

Q1 = ['Do you wish to use; ' 

' 1. your own MATLAB formatted row vector, or ' 

' 2. a program generated vector from the following menu?' 

' - Sine wave ' 

' - Pseudo-noise (PN) sequence ' 

' - Sine wave modulated by PN sequence ' ] ; 
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disp(Ql) 

Pick = input ( 'Answer 1 or 2 : '); 

if Pick == 1, 

% Read in user's input vector 

Wavedata = input ('What is the name of your input vector? '); 
sampfreq = input ('What was the sampling frequency (Hz)? '); 

Tsarople = 1/sampfreq; 
else 

[Wavedata ,Tsample] = wvinput (Pick) ; 
end 

% 

% coefficients are derived from the "h*’ coefficients 

hcoeff = sqrt (2 ) *hcoef f ; 
gcoeff = f 1 iplr (hcoeff ) ; 
for j = 2:2:hlength 

gcoeff (j) = -gcoeff(j); 

end 

% 

% Decomposition algorithm 

datlngth = length (Wavedata) ; 

numrows = ceil ( log (datlngth) /log (2) ) ; % find number of resolution levels 

Lastpts = zeros ( 1 , numrows+1) ; Shift = ones ( 1 , numrows+1 ) ; 

Lastpts (numrows+1 ) = datlngth; lastpt = datlngth; 

for k = numrows : -1 : 1 , % find number of points required at each level 

evnorodd = rem (lastpt , 2 ) ; 
if evnorodd == 0^ 

lastpt = lastpt/2+ceil ( (hlength-2 ) /2 ) ; 
else 

lastpt = ( lastpt-1 ) /2+f ix (hlength/2 ) ; 
end 

Lastpts (k) = lastpt; 

Shift(k) = 2^ (numrows-k+1 ) ; 
end 

hhalf = ceil (hlength/2 ) ; 

dprime = zeros (numrows , Lastpts (numrows) +hhalf) ; 
cprime = zeros (numrows ^ Lastpts (numrows) +hhal f) ; 

clastvct = zeros ( 1 , datlngth+2 *hlength-3 ) ; 
clastvct (hlength-1 : datlngth+hlength-2) = Wavedata; 
ctr = numrows; 
while ctr > 0, 

lastpt = Lastpts (ctr) ; 
nwcvctr = zeros ( 1 ^ lastpt) ; 
nwdvctr = zeros ( 1 , lastpt) ; 

for k = 1: lastpt % coefficients calculated, by resolution level, 

startpt = 2*k-l; % using convolution operation with shifts of 2 

endpt = 2*k+hlength-2 ; % (downsampling) 

nwcvctr(k) = hcoeff *clastvct (startpt : endpt) ' ; 
nwdvctr (k) = gcoef f*clastvct (startpt : endpt) ' ; 

end 

clastvct = [zeros ( 1 , hlength-2) nwcvctr zeros ( 1 , hlength-1 )] ; 
cprime (ctr , 1 : lastpt) = nwcvctr; 
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nwdvctr ; 



dprime (ctr , 1 : lastpt) = 
ctr = ctr-1; 
end 



% 

% Coefficient energies are used as a measure of response 



cenergy 

denergy 

N1 = [ ' 



clc 

disp (b) 
disp (Nl) 
pause 



= cprime.^2; 

= dprime . ^2 ; 

The next plot will show the normalized energy distributions of the' 
approximation and detail coefficients as a function of resolution ' 
level. Based on this information, you will select the minimum ' 

resolution level to be displayed on all future plots. ' 

/ 

'•Edge effects" can cause coefficients to be generated for sample ' 
numbers beyond the endpoints of your data sequence (the lower the ' 
resolution level, the greater the sample number required; therefore, ' 

zero padding of the original data sequence is required) . ' 

/ 

This program translates the indices of the "h" and "g" vectors to ' 
[ -(length of vector-2), 1 ]. As a result, coefficients will exist ' 

for sample numbers beyond the last data point (n=N) , but will not ' 

exist for sample numbers prior to the first data point (n=0) . ' 

These coefficients are necessary for completeness, but may be ' 

difficult to interpret because the number of original data points ' 

used in their calculation can be a small percentage of the total ' 

considered (the rest are "0") . Similarly, coefficients calculated ' 
for sample numbers near the beginning of the sequence will also be ' 
affected by leading zeros affixed to the data sequence. ' 

Limiting the sample numbers required for display may assist inter-' 
pretation; select the lowest resolution level containing significant ' 
energy for this effect. <Return> 



originlE = zeros ( 1 , datlngth) ; originlE = Wavedata.^2; 

Enrgytot = sum (originlE) ; 

Enrgycro = zeros ( 1 , numrows+1 ) ; Enrgydro = zeros (1 ,numrows+l) ; 
ckerr - zeros ( 1 , numrows) ; 

for k = 1: numrows % find energy in each resolution level 

Enrgycro (k) = sum (cenergy (k, :)) ; 

Enrgydro (k) = sum (denergy (k, :)) ; 

end 

if Enrgytot == 0, 

Enrgycro (numrows+1) = 0; 
else 

Enrgycro (numrows+1) = Enrgytot; 

Enrgycro = Enrgycro/Enrgytot ; Enrgydro = Enrgydro/Enrgytot ; % normalize 

end 

Enrgydro (numrows+1) = 0; 



ckenergy = Enrgycro ( 1 ; numrows) +Enrgydro ( 1 : numrows) ; 

% Energy balance check. The total energy at any level should equal the 
% energy in the "c" coefficients in the next higher level, 
for k = 1; numrows 
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ckerr(k) = abs ( Enrgycro (k+1 ) -ckenergy (k) ) ; 
if ckenergy(k) > 10^ (-6) *Enrgytot , 
errenrgy = ckerr (k) /max (Enrgycro (k+ 1) , ckenergy (k) ) ; 
if errenrgy > 10^ (-6), 
clc 

disp(b) 

disp(' ERROR: Energy check between levels does not balance.') 
disp(' If you entered your own "h" vector, recheck its validity.') 

break 
end 
end 
end 

Ivl = [- (numrows) : 1 : 0] ; 

V = [-(numrows+l) , 1 , 0, 1 . 2 ] ; 
axis(v) ; 

subplot(211) , bar ( Ivl , Enrgycro) 

title ( 'Normalized Energy of Approximation Coefficients') 
xlabel ( ' Resolution Level ' ) 
axis (v) ; 

subplot(212) , bar ( Ivl , Enrgydro) 

title ( 'Normalized Energy of Detail Coefficients') 

xlabel ( 'Resolution Level') 

subplot 

pause 

% 

% Output the number of points required to properly display coefficients at 
% each resolution level (determined by "edge effects") 

clc 

disp (b) 

N2 = [ ' Listed below are the number of samples required to properly display' 

' each Resolution Level. ']; 

disp(N2) 
disp (b) 

disp(' Resolution Level Number of samples') 
for k = 1: numrows 

vctrindx = numrows-k+1; 

numsamps = (Lastpts (vctrindx) -1) *Shift (vctrindx) +1 ; 

disp([' ' , num2str (-k) , ' ', num2str (numsamps) ] ) 

end 

nmbrlvl = -input ('What is the lowest Resolution Level you desire to see? '); 

% 

% Plot "c" and "d" coefficients by resolution level 

indxtoO = [ 0 : datlngth-1 ] ; 

plotmin = 1 . 2 *min (Wavedata ) ; plotmax = 1 . 2*max (Wavedata) ; 
if ( (plotmin==0) & (plotmax==0) ) , plotmax = 0.5; end 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 

V = [ 0 , datlngth-1 , plotmin, plotmax] ; 
axis(v) ; 

plot ( indxtoO, Wavedata , '*') 

title ( 'Approximation Coefficients at resolution O') 

xlabel ([ 'Sample number "n" (Time = n * ' ,num2str (Tsample) , ' sec)']) 
pause 

lastrow = numrows-nmbrlvl+l ; 
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clmnsize = (Lastpts ( lastrow) -i) *shif t ( lastrow) +1 ; 

c = zeros (nmbrlvl+1, clmnsize) ; 

c (nmbrlvl+1 , 1 : datlngth) = originlE; 

d = zeros (nmbrlvl , clmnsize) ; 

ctr = numrows; Ctrl = nmbrlvl ; 

while ctr >= lastrow, 

sampval = zeros ( 1, Lastpts (ctr) ) ; 
for k = 1 : Lastpts (ctr) 

index - Shift (ctr )* (k-1) +1 ; 
c (Ctrl , index) = cenergy(ctr,k) ; 
d (Ctrl , index) = denergy(ctr,k) ; 
sampval (k) - index-1; 

end 

plotmin = 1 . 2*min (min (cprime (ctr , : ) ) ,min (dprime (ctr , : ) ) ) ; 
plotmax = 1. 2*max(max(cprime (ctr, : ) ) ,max(dprime(ctr, : ) ) ) ; 
if ( (plotmin==0) & (plotmax==0) ) , plotmax == 0.5; end 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 
V = [0,max(sampval) , plotmin, plotmax] ; 
n = numrows-ctr+1 ; 
axis(v) ; 

subplot (211) , plot (sampval , cprime (ctr, 1: Lastpts (ctr) ),'*') 

title ([ 'Approximation Coefficients at resolution level ' , num2str ( -n) ] ) 

xlabel ([ 'Sample number "n" (Last sample = ', num2str (max (sampval )),')'] ) 

subplot(212) , plot (sampval , dprime (ctr , 1 : Lastpts (ctr) ),'*') 

title ([' Detail Coefficients at resolution level ' , num2str (-n) ] ) 

xlabel ([' Sample number "n" (Time = n * ' , num2str (Tsample) , ' sec)']) 

pause 

clg 

ctr = ctr-1; Ctrl = ctrl-1; 

end 

subplot 



% Reconstruction algorithm 

N3 = [ ' Recomposition of the original data sequence can be checked if ' 

' desired. The recomposition algorithm starts with the approximation' 

.' coefficients at the lowest resolution level previously selected ' 

' and successively "adds** in the detail coefficients of each level ' 

' until the original resolution of the input sequence is obtained. ']; 

clc 

disp(b) 

disp(N3) 

ckalgthm = input ('Do you want to see the Recomposition (Y or N)? ','s'); 
if ckalgthm == 'Y', 

hrecomp = f liplr (hcoef f ) ; grecomp = f liplr (gcoef f ) ; % invert in time 

cstart = cprime (lastrow, :) ; 
ctr = lastrow; 
while ctr <- numrows; 
lastpt = Lastpts (ctr) ; 
ckvctr = zeros(2 , lastpt) ; 

for k = 1: f ix (hlength/2 ) % compute higher level "c" coefficients 

ckvctr(l,:) = ckvctr ( 1 ,:) +hrecomp (2 *k) *cstart (k: lastpt+k-1) +. . 

grecomp (2*k) *dprime (ctr , k: lastpt+k-1) ; 
ckvctr (2 , : ) = ckvctr (2 , : ) +hrecomp (2*k-l) *cstart (k: lastpt+k-1) +. . 

grecomp (2*k-l) *dprime (ctr, k: lastpt+k-1) ; 

end 

if rem(hlength, 2) -= 0, 
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ckvctr(2, : )=ckvctr (2, : ) +hrecoinp (hlength) *cstart (hhalf : lastpt+hhalf-1) + . . 
grecomp (hlength) *dprime (hhalf : lastpt+hhalf-1) ; 

end 

ckvctr = reshape (ckvctr, 1 , 2*Lastpts (ctr) ) ; 
cstart = ckvctr; 

lastsamp = (Lastpts (ctr+1) -1) *Shift (ctr+1) ; 
indxtoO = [ 0 : Shift (ctr+1) : lastsamp] ; 

plotmin = 1 . 2*inin (ckvctr) ; plotmax = 1 . 2*inax (ckvctr ) ; 
if ( (plotmin==0) & (plotmax==0) ) , plotmax = 0.5; end 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 
v = [ 0 , lastsamp, plotmin, plotmax] ; 
axis(v) ; 

plot ( indxtoO , ckvctr ( 1 : Lastpts (ctr+1 ) 

title ([' Level ' , num2str (ctr-numrows) , ' Recomposition']) 

xlabel ([ 'Sample number "n" (Last sample = ', num2str ( lastsamp) ,')'] ) 

pause 

clg 

ctr = ctr+1; 
end 
axis ; 

plot ( indxtoO , ckvctr ( 1 : datlngth) -Wavedata , ' * ' ) 
title (' Recomposition Error') 

xlabel ([' Sample number "n" (Last sample = ', num2str ( lastsamp) ,')'] ) 
pause 
clg 
else 
axis ; 
end 

% 

% Display mesh and contour plots of coefficient energy with optional "zoom" 

% capability to aid analysis 

strtsamp = 1; endsamp = clmnsize; strtrow = 1; endrow = nmbrlvl; zoom = 1; 
highres = -1; lowres = -nmbrlvl; 
while zoom == 1, 

rangea = [ -highres-1 : -lowres ] ; ranged = [ -highres : -lowres] ; 

indxtoO = [strtsamp-1 : endsamp-1] ; 

mesh (c (strtrow: endrow+1 , strtsamp: endsamp) ) 

title (' Energy in Approximation Coefficients') 

if choice - - 1, 

xlabel ([' using Daubechies Wavelet of order ', num2str (choice) ] ) 
end 
pause 

contour (c (strtrow: endrow+l , strtsamp: endsamp) ,10, indxtoO, rangea) 
title (' Contour Map of Approximation Coefficient Energy Distribution') 
xlabel ([' Sample number "n" (Time = n * ' , num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= -Resolution Level)') 
pause 

mesh (d (strtrow: endrow, strtsamp: endsamp) ) 
title (' Energy in Detail Coefficients') 
if choice -= 1, 

xlabel ([ 'using Daubechies Wavelet of order ', num2str (choice) ] ) 
end 
pause 

contour (d (strtrow: endrow, strtsamp: endsamp) , 10, indxtoO, ranged) 
title (' Contour Map of Detail Coefficient Energy Distribution') 
xlabel ([ 'Sample number "n" (Time = n * ', num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= -Resolution Level)') 
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pause 

clc 

showmore = input ('Would you like to "zoom" in on a section (Y or N)? ','s'); 
if showmore == 'Y', 
disp (b) 

disp(' You will set the display sample number and resolution level limits.') 
disp(b) 

disp([' Sample range 0: ' , num2str (clmnsize-1) ] ) 

disp([' Resolution range -1: ' , num2str (-nmbrlvl) ] ) 

revu = input ('Need to see the original contour plots again (Y or N)? ','s'); 
if revu == 'Y', 

contour (c, 10, [ 0 ; clronsize-1] , [0:nmbrlvl] ) 

title ( 'Contour Map of Approximation Coefficient Energy Distribution') 
xlabel ([ 'Sample number "n" (Time = n * ' , num2str (Tsample) , ' sec)']) 
ylabel ( 'Decomposition Number (= -Resolution Level)') 
pause 

contour (d, 10, [ 0 ; clmnsize-1 ] , [l:nmbrlvl] ) 

title ( 'Contour Map of Detail Coefficient Energy Distribution') 
xlabel ([' Sample number "n" (Time = n * ', num2str (Tsample) , ' sec)']) 
ylabel ( 'Decomposition Number (= -Resolution Level)') 
pause 
end 

strtsamp = input ( 'Sample start point? ')+!; 
endsamp = input (' Sample endpoint? ')+!; 

highres = input (' Highest resolution level (least negative)? '); 
lowres = input ( 'Lowest resolution level (most negative)? '); 
endrov = nmbrlvl+1+highres ; 
strtrow = nmbrlvl+1+lowres ; 
else 

zoom = 0 ; 
end 
clg 
end 
clc 
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APPENDIX E 



PROGRAM FOR MULTIPHASE ANALYSIS WITH ANY WAVELET 
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% This program will decompose the input waveform with a Haar wavelet, a 
% Daubechies compactly supported wavelet of user-defined order, or with a 
% wavelet provided by the user. 

% A "sliding” data window is used to maximize signal analysis features 
% of the decomposition by negating the effects of random signal time of 
% arrival. Approximation and detail coefficients may therefore be 
% computed at all resolution levels for every "n" (sample number) . 

% Because of finite input data length, however, some of the coefficients 
% will suffer from "edge effects", i.e., the sample values beyond the 
% end of the data sequence are treated as "0"s. 

% Output graphs consist of "3-D" MATLAB mesh plots and "waterfall" 

% contour maps of the approximation and detail coefficients, in addition 
% to plots of the coefficients at each level (if desired) . 

% 

% Input the "h" coefficients 

clc 

b = [ ' ' ] ; 
disp (b) 

D1 = [ ' Would you like to use: ' 

' 1. your own set of decomposition coefficients, or ' 

' 2. those associated with the Haar wavelet, or ' 

' 3. those derived from Daubechies group of wavelets?']; 

disp(Dl) 

hpick = input ( 'Answer 1, 2, or 3 : '); 
disp (b) 
disp (b) 

if hpick == 1, 

disp(' Note: Sum of coefficients must be normalized to equal 1.') 
disp (b) 

hcoef f=input ( ' What is the name of your decomposition coefficient vector? '); 
choice = 1; 
elseif hpick == 2, 
hcoeff = [0.5 0.5] ; 
choice = 1; 
else 

choice = input ('What is the desired wavelet order (2-10 are available)? '); 
[hcoeff] = daubdata (choice) ; 
end 

hlength = length (hcoef f) ; 

if hlength == 0 
disp(b) 

disp('ERROR: Wavelet order selected is outside allowable range.') 
break 
end 

% 

% Determine if user wants to see plots of the basis functions 

picture=input ( 'Do you want to see the Scaling Function/Wavelet (Y or N)? ','s'); 
if picture == 'Y' 
basisplt 
end 
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% Get the input data vector 



clc 

disp (b) 

Q1 = ['Do you wish to use: ' 

' 1. your own MATLAB formatted row vector, or ' 

' 2. a program generated vector from the following menu?' 

' - Sine wave ' 

' - Pseudo-noise (PN) sequence ' 

' - Sine wave modulated by PN sequence ']; 

disp(Ql) 

Pick = input ( 'Answer 1 or 2 : '); 

if Pick == 1, 

% Read in user's input vector 

Wavedata = input ('What is the name of your input vector? '); 
sampfreq = input ('What was the sampling frequency (Hz)? '); 

Tsample = 1/sampfreq; 
else 

[Wavedata , Tsample] = wvinput (Pick) ; 
end 

I 

% "g" coefficients are derived from the "h" coefficients 

hcoeff = sqrt (2 ) *hcoef f ; 
gcoeff = fliplr (hcoeff ) ; 
for j = 2:2;hlength 

gcoeff (j) = -gcoeff(j); 

end 

% 

% Decomposition algorithm 

datlngth = length (Wavedata) ; 

numrows = ceil ( log (datlngth) /log ( 2) ) ; % determine number of resolution levels 

Lastpts == zeros ( 1 , numrows+1) ; 

for k = 1: numrows % determine number of points required for each level 

Lastpts (numrows-k+l) = datlngth+ (2^ (k) -1) * (hlength-1) ; 

end 

Lastpts (numrows+1) = datlngth; numpts = Lastpts(l); 
d = zeros (numrows , numpts) ; c = zeros (numrows+1 , numpts) ; 
c (numrows+1 , 1 ; datlngth) = Wavedata; cworkvct = Wavedata; 
ctr = numrows; 
while ctr >= 1, 
n = numrows-ctr; 

shift = 2^n; oldendpt = Lastpts (ctr+1) ; newendpt = Lastpts (ctr) ; 
dnewrow = zeros (1 , newendpt) ; cnewrow = zeros ( 1 , newendpt) ; 

for k = 0:hlength-l % coefficient vectors are calculated for each 

% resolution level by adding shifted, weighted 
% versions of the next higher level "c" vector 
% (same effect as direct convolution) 

shiftl = k*shift; 

cnewrow ( 1+shiftl : oldendpt+shiftl) = cnewrow(l+shiftl: oldendpt+shiftl) + . . 

hcoeff (hlength-k) *cworkvct; 

dnewrow ( 1+shiftl : oldendpt+shiftl) = dnewrow(l+shiftl : oldendpt+shiftl) +. . 

gcoeff (hlength-k) *cworkvct ; 

end 

d (ctr , 1 : newendpt) = dnewrow; 
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c (ctr , 1 : newendpt ) = cnewrow; 
cworkvct = cnewrow; 
ctr = ctr-1; 
end 

% 

% User determines output format 

clc 

disp (b) 

D3 = [' Do you want to see: ' 

' 1. separate plots of the coefficients at each Resolution Level' 

' in addition to the "3-D" and contour plots, or ' 

' 2. only the "3-D" and contour plots? ']; 

disp(D3) 

pickplot = input ( 'Answer 1 or 2 : '); 



% Plot "c" and "d" coefficients by resolution level, if desired. 

if pickplot == 1, 

indxtoO = [ 0 : datlngth-1] ; 

plotmin = 1 . 2*min (Wavedata) ; plotmax = 1 . 2*max (Wavedata) ; 
if ( (plotmin==0) & (plotmax==0) ) , plotmax = 0.5; end 
if plotmin > 0, plotmin = 0; end 
if plotmax < 0, plotmax = 0; end 
V = [ 0 , datlngth-1 , plotmin, plotmax] ; 
axis (v) ; 

plot ( indxtoO , Wavedata , ' * ' ) 

title ( 'Approximation Coefficients at resolution O') 

xlabel ([' Sample number "n" (Time = n * ' , num2str (Tsample) , ' sec)']) 
pause 

for k = numrows:-l:l 
n = numrows-k+1; 
endpt = Lastpts(k); 
indxtoO = [0:endpt-l]; 

plotmin = 1 .2*min(min(c(k, : ) ) ,min(d(k, : ) ) ) ; 

plotmax = 1 . 2*max (max (c (k, : ) ) , max (d (k, : ) ) ) ; 

if ( (plotmin==0) & (plotmax==0) ) , plotmax = 0.5; end 

if plotmin > 0, plotmin = 0; end 

if plotmax < 0, plotmax = 0; end 

V = [0,endpt-l , plotmin, plotmax ] ; 

axis (V) ; 

subplot (211) , plot (indxtoO , c (k, 1 : endpt) , ' * ' ) 

title ([ 'Approximation Coefficients at resolution level ' , num2str (-n) ] ) 
xlabel ([' Sample number ”n" (Number of points = ', num2str (endpt ),')'] ) 
subplot (212 ) , plot (indxtoO , d (k, 1 : endpt) , ' * ' ) 

title ([' Detail Coefficients at resolution level ' , num2str (-n) ] ) 
xlabel ([ 'Sample number "n" (Time = n * ', num2str (Tsample) , ' sec)']) 
pause 
clg 
end 

subplot 
axis ; 
end 



% Coefficient energies are used as a measure of response. Display mesh and 
% contour plots of coefficient energy with optimal "zoom*' capability to aid 
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% signal analysis 



c = c . ^ 2 ; 
d = d.-2; 

strtsamp = 1; endsainp = numpts; strtrow = 1; endrow = numrows; zoom = 1; 
highres = ~1; lowres = -numrows; 
while zoom == 1, 

rangea = [ -highres-1 : -lowres] ; ranged = [ -highres: -lowres] ; 

indxtoO = [ strtsamp-1 : endsamp-1 ] ; 

mesh (c (strtrow: endrow+l, strtsamp: endsamp) ) 

title ( 'Energy in Approximation Coefficients') 

if choice -= 1, 

xlabel ([ 'using Daubechies Wavelet of order num2str (choice) ] ) 
end 
pause 

contour (c( strtrow: endrow+l, strtsamp : endsamp) , 10, indxtoO , rangea) 
title (' Contour Map of Approximation Coefficient Energy Distribution') 
xlabel ([ 'Sample number "n" (Time = n * ' , num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= - Resolution Level)') 
pause 

mesh (d (strtrow: endrow, strtsamp : endsamp) ) 
title ( 'Energy in Detail Coefficients') 
if choice -= 1, 

xlabel ([' using Daubechies Wavelet of order ', num2str (choice) ] ) 

end 

pause 

contour (d (strtrow: endrow, strtsamp: endsamp) , 10, indxtoO, ranged) 

title (' Contour Map of Detail Coefficient Energy Distribution') 

xlabel ([' Sample number "n” (Time = n * ' , num2str (Tsample) , ' sec)']) 

ylabel (' Decomposition Number (= - Resolution Level)') 

pause 

clc 

showmore = input ('Would you like to "zoom” in on a section (Y or N)? ','s'); 
if shovTHore == 'Y' , 
disp (b) 

disp(' You will set the display sample number and resolution level limits.') 
disp (b) 

disp([' Sample range 0 : ' , num2str (numpts-1) ] ) 

disp([' Resolution range -1 :', num2str (-numrows) ] ) 

revu = input ('Need to see the original contour plots again (Y or N)? ','s'); 
if revu == ' Y ' , 

contour (c, 10, [0:numpts-l] , [ 0: numrows] ) 

title (' Contour Map of Approximation Coefficient Energy Distribution') 
xlabel ([' Sample number "n" (Time = n * ', num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= - Resolution Level)') 
pause 

contour (d, 10, [ 0:numpts-l] , [ l:numrows] ) 

title (' Contour Map of Detail Coefficient Energy Distribution') 
xlabel ([' Sample number "n" (Time = n * ', num2str (Tsample) , ' sec)']) 
ylabel (' Decomposition Number (= - Resolution Level)') 
pause 
end 

strtsamp = input (' Sample start point? ')+l; 
endsamp = input (' Sample endpoint? ')+l; 

highres = input ( 'Highest resolution level (least negative)? '); 
lowres = input (' Lowest resolution level (most negative)? '); 
endrow = numrows+l+highres ; 
strtrow = numrows+l+lowres ; 
else 
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zoom = 0 ; 
end 
clg 
end 
clc 
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APPENDIX F 



PROGRAM FOR GENERATING ITERATIVE PLOTS OF WAVELETS 
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% This program creates and plots the desired iterative approximation 
% to the scaling function and the associated wavelet as determined by 
% the input "h*' vector. The construction is based on the "graphical” 
% recursion method. 



hcoefnew = 2*hcoeff; 
m = length (hcoefnew); 



% The "g" vector is determined from the "h" vector. 

gcoefnew = fliplr (hcoefnew) ; 
for j = 2:2:m 

gcoefnew(j) = -gcoef new ( j ) ; 
end 

% 

% Recursively build basis functions 

numbrits = input ( 'How many iterations shall we run for the approximation? '); 
size = 1; 

hpast = ones (l,m) ; 
newsize = m; 

for i = 1: numbrits 

% Build scaling function approximation 

hnew = zeros ( 1 , newsize) ; 
for k = 0:size-l 

hnew(2*k+l : 2*k+m) =hnew ( 2 *k+l : 2*k+m) +hpast (k+1) *hcoefnew; 

end 

% Get wavelet from scaling function approximation and also build time vector 

wv = zeros (1 , 2*newsize) ;wvlet = zeros (1 , newsize) ;timevctr = zeros ( 1 , newsize) ; 
shift = newsize/ (ro-1) ; 
for k = 0:m-l 

shiftl = round(k*shift+l) ; 
shift2 = round (newsize+k*shift) ; 

wv (shiftl : shift2) = wv(shiftl:shift2)+gcoefnew(k+l) *hnew; 

end 

for k = 1: newsize 

wvlet(k) = wv(2*k-l) ; 
timevctr(k) = k-1; 

end 

interval = (1/2) ^i; 

timevctr = interval*timevctr ; % scale time axis 

s = interval*newsize ; 

% Plot basis functions and check calculations with areas and inner products 
if i <= 5, 

plot (timevctr , hnew, ' , timevctr , wvlet , ' * ' ) 

xlabel ([ 'Support = ' , num2str (s) , ' (Scaling Function: ++++, Wavelet: ****)']) 
else 
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plot ( timevctr , hnew, ' - ' , timevctr , wvlet , ' — ' ) 

xlabel ([' Support = ' , num2str ( s) , ' (Scaling Function: line, Wavelet: dash)']) 
end 

title ([ 'Approximations for Iteration number ' , num2str ( i ) ] ) 
pause 

hpast = hnew; 
size = newsize; 
newsize = 2*size-^m-2; 



phisum = sum (hpast) ; 
wvsum = sum(wvlet) ; 
sp(i) = phisum*interval ; 
sw(i) = wvsum* interval ; 
ip(i)= hpast*wvlet' ; 



% find area under scaling function 
% find area under wavelet 

% find scaling function/wavelet inner product 



end 



% 

% Output calculation checks 



disp(' ') 

disp(' Area under Area under 

disp ( 'Scaling Function Wavelet 
disp(' ') 

for i = l:numbrits 

disp([' ' ,num2str (sp(i) ) , ' 

end 

disp(' ') 

disp ( ' <Return> ' ) 

pause 



Inner ' ) 
Product ' ) 



' , num2str (sw ( i) ) 



clear hcoefnew gcoefnew hpast timevctr w 



,num2str (ip(i) ) ]) 
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APPENDIX G 



FUNCTION FOR CONSTRUCTION OF INPUT DATA 
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function [ Data , Tsample ] = wvinput(x) 

% This function supports the various multiresolution programs by building 
% the input test vector desired by the user from the following choices: 

% Sine wave, Pseudo-noise sequence, and Sine wave modulated by Pseudo- 
% noise sequence w^ith optional additive Gaussian noise. 



% Determine desired input signal type, 
b = [ ' ' ] ; 

Q2 = [ ' Select the desired type of waveform: ' 

' 1 . Sine wave ' 

' 2. Pseudo-noise (PN) sequence ' 

' 3. Sine wave modulated by PN sequence']; 

disp (b) 
disp(Q2) 

Make = input ( 'Answer 1, 2 or 3 : '); 



% Generate the desired sinusoidal signal, 



if Make -= 2, 
clc 

% Determine desired sine wave characteristics. 

N2 = [ ' You will determine the following Sine wave parameters:' 

' - frequency ' 

' - phase ' 

' - amplitude ' 

' - number of samples ' 

' - sampling frequency or data length ']; 

disp (b) 
disp(N2) 

fc = input ( 'Sinewave frequency ("fc") in Hertz? '); 
theta = input ('Phase (degrees)? '); 

A = input ( 'Amplitude? '); 

N = input ( 'Number of Samples ("N”, with N = a power of 2)? '); 

Q3 = [' Do you wish to set: ' 

' 1. the sampling frequency, or' 

' 2. the data length? ']; 

disp (b) 
disp(Q3) 

S5 = input ( 'Answer 1 or 2 . '); 



if S5 == 1, 

N3 = [ ' Note: Sampling frequency (fs) must be selected so that ' 

' fs/fc >= 2; Data length will be set to allow "N" samples.']; 
disp (b) 
disp (N3 ) 

fs2fc = input ( 'Desired sampling-to-signal frequency ratio? '); 
const = 2*pi/fs2fc; 
fs = fs2fc*fc; 

else 

N4 = [ ' Note: Sampling frequency will be chosen to give ”N*‘ samples' 

' so observation time must be <= "N"/ ( 2 * ” f c" ) . ']; 

disp (b) 
disp(N4 ) 
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Tobserv = input (' Desired observation time (seconds)? '); 
fs = N/Tobserv; 
const = 2*pi*fc/fs; 

Tsample = 1/fs; 

end 

Tsample = 1/fs; 

% Construct sine wave 
theta = theta*pi/180 ; 
for k = 1:N 

Data(k) = sin (const* (k-1) +theta) ; 

end 

Data = A* Data ; 

if Make == 3 , % used if PN modulation is desired 

Sinedata - Data; 
disp (b) 
disp (b) 

N4pt5 = [' You will now determine the PN sequence characteristics.']; 
disp(N4pt5) 

end 

end 



% This section generates the desired pseudo-noise signal by building the 
% appropriate feedback shift register (FSR) . 

if Make -= 1, 
clc 

% Determine desired PN sequence characteristics. 

N5 = [ ' Note: PN sequence = (+1 or -l,...,2^m - 1} ' 

' Feedback shift register initial state is { 1 , 1 , . . . , 1 ) ' ] ; 
disp (b) 
disp(N5) 

N6 = [ ' You will determine the following parameters:' 

' - number of stages ' 

' - number and position of taps ' 

' - chip rate ' 

' - time delay ' 

' - number of samples ' 

' - sampling frequency or data length ']; 

disp (b) 
disp(N6) 

mstages = input (' Desired number (m) of stages (2-10)? '); 
numbrtap = input ( 'Number of taps? '); 
for k = 1: numbrtap 

Tap(k) -input ([ 'What is the position of Tap number ' , num2str (k) , ' ? ']); 

end 

chiprate = input('Chip rate (chips/sec)? '); 

delay = input (' Sequence delay in chips (positive real number)? '); 
if Make == 2, 

% If the PN signal alone is desired the user must input the number 
% of samples and sampling rate; otherwise, they are dictated by 
% the choices made for the sinusoid. 

A = sqrt (2) ; 

N = input (' Number of samples? '); 

Q4 = [' Do you wish to set: ' 
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' 1. the sampling frequency, or' 

' 2. the data length? ']; 

disp (b) 
disp(Q4) 

P6 = input ( 'Answer 1 or 2 . '); 
if P6 “ 1, 

N7 = [ ' Note: Sampling frequency ("fs") must be selected so that ' 

' fs/chiprate >= 2; data length will be set for *'N" samples.'] 
disp (b) 
disp(N7) 

fs2cr = input (' Desired sampling frequency to chip rate ratio? '); 
Tsample = 1/ (chiprate*f s2cr) ; 
else 

N8 = [ ' Note: Sampling frequency will be chosen to give **N" samples' 
' so observation time must be <= *'N**/ (2*Chip rate). ' 

disp (b) 
disp(N8) 

Tobserv = input (' Desired observation time (seconds)? '); 
fs = N/Tobserv; 
fs2cr = fs/chiprate; 

Tsample = 1/fs; 
end 
else 

disp(b) 

N9 = [ ' Same number of samples ( ' , num2str (I^) , ' ) and sampling 
NIO = [' frequency ( ' , num2str ( fs) , ' Hz) that were selected for']; 

Nil = [' the sine wave will be used. <Return>']; 

disp (N9) 

disp(NlO) 

disp(Nll) 

pause 

fs2cr - fs/chiprate; 
end 



% Generate base PN sequence 
FSR = ones ( 1 , mstages) ; 

L = 2 ^mstages - 1; 
cklength = N/(fs2cr*L); 

repeat = ceil (cklength) ; % Determine how many periods are necessary 

V = zeros (1 , (repeat+1) *L) ; 

for k = 1:L % build one period of the sequence 

V(k) = FSR (mstages) ; 
sigma = 0; 

for m = linumbrtap % modulo-2 tap adder 

sigma = sigma+FSR (Tap (m) ) ; 

end 

FSR(2 :mstages) = FSR ( 1 :mstages-l) ; 

FSR(l) = rem (sigma , 2 ) ; 

end 

% Ensure enough periods are available for the desired number of samples 
for n = 1: repeat 

V(n*L+l: (n+1) *L) =V(1:L); 

end 

delay = rem (delay , L) ; 
for n = 1:N 

Data(n) = V ( f ix ( (n-1) /fs2cr + delay)+l); 
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if Data(n) == 0, Data(n) 



-1; end 



% make sequence bipolar 



end 

end 

% 

% If selected, modulate the sine wave with the PN sequence. 

if Make == 3 

Data = Sinedata . *Data ; 
end 

% 

% Plot the constructed signal and add noise, if desired. 
indxtoO = [0:N-1]; 

V = [ 0 , N-1, -1 . 2*max (Data) , 1 . 2*max (Data) ] ; 
axis(v) ; 

plot (indxtoO, Data) , title ( 'MATLAB Plot of Input Signal Vector') 
xlabel (' Sample number') ,ylabel ( 'Volts' ) 
pause 
axis ; 

clc 

disp (b) 

disp(' White Gaussian Noise is generated by the MATLAB "random" function.') 
noise = input ('Do you want noise added to your signal (Y or N) ? ','s'); 
if noise == ' Y ' , 

% SNR based on "continous" signal energy and the Gaussian noise variance 

SNRdB = input('What is the desired Signal-to-Noise Ratio (dB)? '); 

SNR = 10^ (SNRdB/10) ; 

Noisepwr = ( (A^2 ) /2 ) /SNR ; 
rand (' seed ', 0) ; rand ( ' noraal ' ) ; 
noisvctr = Noisepvr*rand ( 1 , N) ; 

Data = Data+noisvctr ; 

V = [ 0 , N-1 , -1 . 2*max (Data) , 1 . 2*max (Data) ] ; 
axis(v) ; 

plot (indxtoO , Data) , title ( 'MATLAB Plot of Input Data Vector') 
xlabel ( 'Sample number' ) ,ylabel ( 'Volts' ) 
pause 
axis ; 
end 



return 
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APPENDIX H 



FUNCTION FOR DAUBECHIES* H-COEFFICIENTS 
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function [hcoeff] = daubdata(x) 

% This function returns the proper "h" coefficients for the specified order 
% Daubechies wavelet. The input value of *'x" is the specified order. 



if X == 2 
hcoeff = 
elseif X == 
hcoeff = 

elseif X == 
hcoeff = 

elseif X == 
hcoeff = 



elseif X == 



[ 0.4 829 629 13 14 5; 0.83 6516303 73 8; 0.224 14 38 68 04 2 ;-0. 129409522551] ; 

3 

[ 0.332 67 05529 50 JO. 8 0689 15093 11; 0.4 59877 502118; -0.135011020010 ; 
-0. 085441273882; 0.035226291882] ; 

4 

[0.230377813309 ; 0 . 7 1484 657 0553 ; 0 . 63 08 807 679 3 0 ; -0 . 02798 3 7 6 94 17 ; 
-0. 187034811719 ;0. 030841381836 ;0. 032883011667 ; -0 . 0105974 01785 ] ; 

5 

[0.160102397974 ; 0 . 603829269797 ;0.724308528438;0.138428145901; 

-0. 24 2294 887066; -0.03 224 4 869 585; 0. 077 5714 93 840; -0. 006241490213 ; 
-0. 012580751999 ;0. 003335725285] ; 

6 



hcoeff 



elseif X 
hcoeff 



elseif X 
hcoeff 



elseif X 
hcoeff 



elseif X 
hcoeff 



else 



= [0.111540743350;0.494623890398;0.751133908021;0.315250351709; 

-0. 22 62 64 69 3965; -0. 129766867567 ; 0 . 097501605587 ;0.027522865530; 
-0. 03 158203 9318; 0. 000553 84 2201; 0. 004777257 511; -0.001077 3 01085] ; 

= [0. 077852054085 ;0. 396539319482 ;0.729132090846;0.469782287405; 
-0.143906003929 ; -0 . 224 03 6184 994 ;0.071309219267;0.080612609151; 
-0. 03 8 02993 693 5; -0. 016574 541631 ;0. 012550998 556; 0. 0004 2957 7973 ; 
-0.001801640704 ; 0. 000353713800] ; 

== 8 

= [0.054415842243 ; 0 . 3 12871590914 ; 0 . 67563 073 6297 ; 0 . 58 53 54 6 8 3 654 ; 

-0. 015829105256 ; -0.284 01554 2962; 0. 000472484574 ; 0 . 12874 74 26620 ; 
-0. 017369301002 ; -0 . 044 08825393 1 ; 0 . 013981027917 ; 0 . 00874 6094047 ; 
-0. 004 87 03 529 9 3 ; -0 . 0003 917 4 0373 ; 0 . 00067 54 4 94 06 ; -0. 000117 47 678 4 ] ; 
== 9 

= [0.038077947364; 0.243834674613 ;0.604823123690;0.657288078051; 
0.133 19738 58 25; -0.293273783279; -0.09 684 07 83223; 0.14 854 07 49338; 
0.030725681479 ; -0 . 06763 2829061 ; 0 . 0002 5 0947 1 15 ; 0 . 022 3 616 62 124 ; 
-0.004723204758 ; -0 . 004281503682 ; 0. 001847646883 ; 0 . 0002 3 03 8 5764 ; 
-0.000251963189 ; 0 . 00003 9 3 4 7 3 2 0 ] ; 

== 10 

= [ 0.02667 0057901 ;0. 188 17 68 0007 8; 0.5272 0118893 2; 0.6884 59 03 9454 ; 

0. 28 117234 3 661; -0.2 4 984 64 24 3 27 ; -0 . 19594 627 4 3 77 ; 0 . 1 27 3 69 3 4 03 3 6 ; 
0.093057364604 ; -0 . 0713 94 147 166 ; -0 . 0294 5753 6822 ; 0 . 03 3 2 12 674 059 ; 
0.003606553567 ; -0 . 01073 3 1754 83 ;0.001395351747;0.001992405295; 

-0. 000685856695;-0.000116466855;0. 000093 58 8670; -0. 000013264203] ; 



hcoeff = [ ] ; 
break 
end 



hcoeff = hcoef f Vsqrt (2) ; % normalize the "h" vector 

return 
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